Amplification and detection of single molecule conformational fluctuation through a 
protein interaction network with bimodal distributions 
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A protein undergoes conformational dynamics with multiple time scales, which results in fluctu- 
ating enzyme activities. Recent studies in single molecule enzymology have observe this "age-old" 
dynamic disorder phenomenon directly. However, the single molecule technique has its limitation. 
To be able to observe this molecular effect with real biochemical functions in situ, we propose 
to couple the fluctuations in enzymatic activity to noise propagations in small protein interaction 
networks such as zeroth order ultra-sensitive phosphorylation-dephosphorylation cycle. We showed 
that enzyme fluctuations could indeed be amplified by orders of magnitude into fluctuations in 
the level of substrate phosphorylation — a quantity widely interested in cellular biology. Enzyme 
conformational fluctuations sufficiently slower than the catalytic reaction turn over rate result in a 
bimodal concentration distribution of the phosphorylated substrate. In return, this network ampli- 
fied single enzyme fiuctuation can be used as a novel biochemical "reporter" for measuring single 
enzyme conformational fiuctuation rates. 
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I. INTRODUCTION 



The dynamics of an enzyme is usually characterized by 
a rate constant describing its catalytic capacity, which is 
a standard practice on studying dynamics of enzymes 
and enzyme-involved networks. Recent advances in sin- 
gle molecule techniques allow examining enzyme activ- 
ities at single molecule levels [1 U, i, i, H, U, S B i, 
E, EH, m, [II- It is found that the rate "constant" 
of an enzyme is in general a broad distribution. Physi- 
cally, it is because the enzyme conformation is under con- 
stant fluctuation at varying time scales [H, [13] . Single 
molecule techniques can measure the instant rate con- 
stants at a given conformation. The single molecule re- 
sults are consistent with extensive early biochemistry and 
biophysics studies. Biochemists have long noticed that 
protein conformational fluctuations (which can be in the 
time scale from subsecond to minutes and even hours) can 
be comparable and even slower than the corresponding 
chemical reactions (usually in the range of subsecond) 
[l5| . Slow conformational motions result in hysteretic 
response of enzymes to concentration changes of regula- 
tory molecules, and cooperative dependence on substrate 
concentrations [l^, [l^, [13, In physical chemistry, 

the term dynamic disorder is used for the phenomenon 
that the rate of a process may be stochastically time- 
dependent p^ . Extensive experimental and theoretical 
studies exist since the pioneering work of Frauenfelder 
and coworkers poj . AUosteric enzymes can be viewed 
as another class of examples. According to the classi- 
cal Monod-Wyman-Changeux model and recent popula- 
tion shift model, an allosteric enzyme coexists in more 
than one conformation [U, . Recent experiments also 
show that the conformational transition of an allosteric 
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enzyme happens in micro- to millisecond time scale or 
longer Xing proposed that in general internal con- 
formational change should be considered on describing 
enzymatic reactions, and it may have possible implication 
on allosteric regulation mechanism 24|. Wei et. al also 
suggested a similar formalism for describing enzymatic 
reactions [1^ . Current single molecule enzymology stud- 
ies focus on metabolic enzymes. It remains an important 
unanswered question if dynamic disorder is a general phe- 
nomenon for enzymes, e.g., the enzymes involved in sig- 
nal transduction. While the technique reveals important 
information, the single molecule approach also has strict 
requirement on the system to achieve single molecule sen- 
sitivity, which limits its usage for quick and large scale 
scanning of enzymes. 

In this work we discuss an idea of coupling molecu- 
lar enzymatic conformational fluctuations to the dynam- 
ics of small protein interaction networks. Specifically 
we will examine a phosphorylation/dephosphorylation 
cycle (PdPC) ^2^1]. Our analysis will be applicable to 
other mathematically equivalent systems, such as GTP- 
associated cycle, or more general a system involving two 
enzymes/enzyme complexes with opposing functions on 
a substrate. As an example for the latter, the system 
can be an enzymatic reaction consuming ATP hydrolysis 
{e.g., a protein motor) coupled to a ATP regeneration 
system-in this case ATP is the substrate. The PdPC is 
a basic functional module for a wide variety of cellular 
communications and control processes. The substrate 
molecules can exist in the phosphorylated and dephos- 
phorylated form, which are catalyzed by kinase and phos- 
phatase respectively at the expense of ATP hydrolysis. 
The percentage of the phophorylated substrate form de- 
pends on the ratio of kinase and phosphatase activities in 
a switch like manner called ultra-sensitivity. Through the 
PdPC, slow conformational (and thus enzymatic activ- 
ity) fluctuations at the single molecule level can be ampli- 
fied to fluctuations of substrate phosphorylation forms by 
several orders of magnitude, and make it easier to detect. 
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The coupling between molecular fluctuations and net- 
work fluctuations itself is an interesting biological prob- 
lem. Recent studies revealed that the intrinsic/extrinsic 
noise, when it is introduced into the biological system, 
has significant influences on the behavior and sensitivity 
of the entire network [23,|2i[2i,[33,l3ll[3i,[3l|33. Un- 
like the noise sources studies previously, the internal noise 
due to dynamic disorder shows broad time scale distri- 
butions. Its effect on network-level dynamics is not well 
studied [s^ . We will show that a bimodal distribution of 
the PdPC substrate form can arise due to dynamic dis- 
order, which may have profound biological consequences. 



II. THE MODEL 

A PdPC is shown in [T^. X and X* are the unphos- 
phorylated and phosphorylated forms of the substrate, 
respectively. We assume A is the phosphatase obeying 
normal Michaelis-Menten kinetics. The kinase E, on the 
other hand, can assume two conformations with differ- 
ent catalytic capacity. In general an enzyme can assume 
many different conformations. The two state model here 
can be viewed as coarse-graining. The set of reactions 
describing the system dynamics are listed below: 

X -\- El ^ X El ^ X -\- El , 

kit fc+b 

X -\' E2 ^ XE2 ^ X + E2, 

^Ib' ^' + b' 

k2 f k— f 

X* +aJ x*a ^ x + a, 

k2b k—h 
a a. 
13 13' 

To ensure proper detailed balance constraint, each pair 
of forward and backward reaction constants are related 
by the relation A^" = -kBTln{kf /h), with AfjP the 
standard chemical potential difference between the prod- 
uct (s) and the reactant(s), the Boltzmann constant, 
T the temperature. Exceptions are the three chemical 
reaction steps, which we assume couple to ATP hydroly- 
sis, and thus extra terms related to ATP hydrolysis free 
energy AfiATP are added, we assume A^°-f = 
—kBT\n{kf/kb) for each of these reactions, so one ATP 
molecule is consumed after one cycle. In general the con- 
former conversion rates a, /3, a' , (3' are different. For sim- 
plicity in this work, we choose a ^ f3 — a' — l3' unless 
specified otherwise. 

Let's define the response curve of the system to be the 
steady-state percentage of X* as a function of the cat- 
alytic reactivity ratio between the kinase and the phos- 
phatase {9 = {pik++p2k'^)NEj{k-NAt), where Pi is the 
probability for the kinase to be in conformer i) , and Ne^ 
and Njit are total numbers of kinase and phosphatase 
molecules. At bulk concentration, for both fast and slow 
kinase conformational change, the response curve shows 



usual sigmoidal but monotonic dependence (see ^p) [2g| . 
Here we only consider the zero-th order regime where the 
total substrate concentration is much higher than that of 
the enzymes. 

However, the situation is different for a system with 
small number of molecules. For simplicity let's focus on 
the case with one kinase molecule. Physically, suppose 
that the average reactivity ratio 6^1 (not necessarily ex- 
actly at 1), but the corresponding 9i = k^/ {Nji^k^) > 1 
and 02 = k'j^/ {N A^k^) < 1. Consequently, the substrate 
conversion reactions are subject to fluctuating enzyme 
activities, an manifestation of the molecular level dy- 
namic disorder. Because of the ultra-sensitive nature 
of the PdPC, small enzyme activity fluctuation (in the 
vicinity of ^ = 1) can be amplified into large fluctuations 
of substrate forms (in the branches of high or low num- 
bers of X*). The relevant time scales in the system are 
the average dwelling time of the kinase at the new confor- 
mation tk , and the average time required for the system 
to relax to a new steady-state substrate distribution once 
the kinase switches its conformation Tg. The former is re- 
lated to the conformer conversion rates. The latter is de- 
termined by the enzymatic reaction dynamics as well as 
the number of substrate molecules. If the kinase confor- 
mational switch is sufficiently slow (tk >> T5), so that 
on average for the time the kinase dwelling on each con- 
formation, the substrate can establish the steady state 
corresponding to 9i, which is peaked at either high or 
low Nx*- Then the overall steady-state substrate dis- 
tribution is a bimodal distribution, which is roughly a 
direct sum of these two single peaked distributions. This 
situation resembles the static limit of molecular disor- 
ders [Tgj . Increasing conformer switching rates tends to 
accumulate population between the two peaks, and even- 
tually results in a single-peaked distribution {tk < ts). 
A critical value of a (or /3) exists where tk ^ T5, and 
one peak of the distribution disappears. There arc two 
sets of {tk, ts) corresponding to the transition from con- 
former 1 to 2, and vice versa. In principle, in the slow 
enzyme conform conversion regime where the substrate 
shows non-unimodal distribution {tk > ts), one can ex- 
tract molecular information of the enzyme fluctuations 
from the greatly amplifled substrate fluctuations. This is 
the basic idea of this work. 



III. NUMERICAL STUDIES 

To test the idea, we performed stochastic simulations 
with the Gillespie algorithm ^36i] using the parameters 
listed in Table 1 and Appendix A. [TJ: shows that with 
a single slowly converting two-state kinase, the number 
of X* jumps between high and low values, and shows bi- 
modal distribution (see [TJi) . [5] gives systematic studies 
on this phenomenon. There exist two critical values of a, 
ai,OL2. With a < (ai, 02), the substrate distribution has 
two well-separated peaks (H^). On increasing a, the two 
original peaks diminish gradually while the the region 
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between the two peaks accumulates population to form 
a new peak ( ^jp) . The two original peaks disappear at 
a = ai and a — a2 respectively, and eventually the dis- 
tribution becomes single-peaked ( [^-e) . [5f summarizes 
the above process using the distance between peaks. The 
results divide into three regions. The point of transition 
between the left and the middle regions indicates disap- 
pearance of the left peak (corresponding to Ht) . That 
between the middle and the right regions indicates dis- 
appearance of the right peak (corresponding to Eli). 

[2^ shows that the critical value of a decreases with 
the total substrate number Nxt- An increased num- 
ber of Nxt gives a larger T5, which requires a larger tk 
(slower conformation conversion rate) in order to gener- 
ate a multi-peaked distribution. [2^ also compares the- 
oretical (see below) and simulated critical a values at 
different values of Nxt ■ The plot shows that the simula- 
tion results agree reasonably with the theoretical predic- 
tions, although the simulated critical a is slightly smaller 
than the theoretical values, which means that the peak 
disappears earlier on increasing fluctuation rate a. The 
discrepancy of the two could be due to stochasticity of en- 
zymatic reactions, which is fully accounted for in the sim- 
ulations, but neglected in the theoretical treatment. The 
broader distribution leads to an earlier disappearance of 
peaks. This argument is supported by the fact that the 
difference between theoretical and simulated critical a is 
getting closer when the substrate number becomes larger, 
so fluctuations due to enzymatic reactions are further 
suppressed. 

In the above discussions, we focus on a system with a 
single copy of the two-state kinase molecule. Appendix 
B shows that a simulation result with multi-state model 
gives similar behaviors. [3^ shows that with multiple 
copies of enzymes, the substrate distribution of a PdPC 
can show similar transition from bimodal (or multi-modal 
in some cases) to unimodal behaviors, but the critical 
values of a are smaller (corresponding to slower confor- 
mational change) than those for the single kinase case . 
In these calculations, we scale the system proportionally 
to keep all the concentrations constant. 

Possible biological significances of the bimodal distri- 
butions will be discussed below. Here we propose that 
additionally one can use the phenomenon to extract sin- 
gle molecule fluctuation information, especially the con- 
former conversion rates. Conventionally the information 
is obtained through single molecule experiments [1, [ll] . 
For simplicity here we focus on the single enzyme case 
only. Suppose that an enzyme fluctuates slowly between 
different conformers, and one can couple a single molecule 
enzyme with a PdPC (or a similar system) with fast en- 
zymatic kinetics. Then the conformational fluctuation 
dynamics at the single molecule level will be amplified to 
the substrate form fluctuations by orders of magnitude. 
|4] gives the result of such an experiment simulated by 
computer. The trajectory clearly shows two states. To 
estimate the time the system dwelling on each state, we 
define the starting and ending dwelling time as the first 



time the number of substrate molecules in the X* form 
reaches the peak value of Nx* distribution corresponding 
to that state in the forward and backward direction of the 
trajectory (see|3^). The above algorithm of finding the 
dwelling time may miss those with very short dwelling 
time so the substrate may not have enough time to reach 
the peak value, as seen in the trajectory. Nevertheless, 
the obtained dwelling time distributions are well fitted 
by exponential functions. The exponents give the values 
of a and /3, in this case ~ 0.8 x 10~^ for both of them, 
which are good estimations of the true value 10^'^. 

IV. THEORETICAL ANALYSIS 
A. Analytical estimate of the critical points 

Here we provide quantitative analysis of the above time 
scale argument. Let's define 

k+[Et]x , k_[At]{[X]t - X 
Km+ + X Km- + [Xt\ - X 

and similar expression for the case that the kinase as- 
sumes conformer 2 except fc+ is replaced by fc'^ . Then the 
kinetics of a PdPC with one two-state kinase molecule 
is governed by a set of Liouville equations under the 
Langevin dynamics approximation (lol . IstI . [ssj , 

dfPiix) = ^9^2; [{gi + hi)pi{x)] 

[{-gi + hi)pi{x)] - api{x) + f3p2{xl^) 

dtP2{x) = ■^d^^[{g2 + h2)p2{x)] 

-dx [(-32 + h2)p2{x)] + api{x) - /3p2(a;](3) 

where pi{x) is the probability density to find the system 
at kinase confomer i and the number of substrate form 
X being x, is the system volume. For mathematical 
simplicity in the following derivations, we assume that 
for a given kinase conformation, the substrate dynam- 
ics can be described continuously and deterministically. 
This approximation is partially justified by the relative 
large number of substrates. Then one can drop the diffu- 
sion term containing Q, and solve the above equations an- 
alytically (see Appendix C). The theoretical steady state 
solutions of p{x) = pi{x) + P2{x) are also plotted in E) 
Since in our analysis we neglected stochasticity of enzyme 
reactions due to the finite number of substrates, the an- 
alytical solutions are bound by the two roots {xi,X2) of 
the equations, 

/i(a;i) =0,/2(x2) = 0. (4) 

In the case of fast switching rates, the solutions vanish at 
the turning points {xi,X2), and become identically zeros 
outside of the interval [xi, X2\- Physically it means that 
the enzyme equilibrates quickly and the rest of the system 
'feels' only averaged reactivity of two enzyme's states. 
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In the regime of slow switching rates, the steady state 
solutions pi{x) and pi{x) have two integratable singu- 
lar points at {xi,X2)- The solutions diverge at these 
points (although integration of pi over x is still finite). 
Of course, the neglected 'diffusion' term due to sub- 
strate fluctuations becomes important in this situation. 
This term smears the singularities at the turning points 
{xi,X2)- Kepler et al. noticed similar behaviors in their 
simulation results. However, in order to estimate the 
critical values of the switching rates ac and f3c, which 
correspond to transition between unimodal and bimodal 
distribution, the set of equations Eq. ^ and Eq. ^ are 
sufficient. 

Even within this approximation the analytical predic- 
tion of the transition points agrees well with the simu- 
lation results. The conditions for disappearance of these 
two single points are, 

a > ac ^ dxfi{xi),P > Pc = dxf2{x2)- (5) 

Note that 1/a (or 1//3) is the average dwelling time of 
an enzyme configuration r^^, and 1/ac (or 1/Pc) is the 
relaxation time after the system linearly deviates from 
the single conformer steady state (at xi or X2), which 
are T5 in the previous discussions. ET & g show good 
agreement between the critical points obtained by sim- 
ulation and by Eq. ([5]). The agreement becomes better 
for larger number of substrates, suggesting that the dis- 
crepancy between the simulated and theoretical results 
are due to neglecting substrate fluctuations in the theo- 
retical treatment. 



Parameters 


Values 


rate constant 


in reduced unit 


kif 


50 


k+f 


1.6 


kif 


50 


k+f, 


0.4 


k2f 


50 


k-f 


1 


Free energy 


in ksT 


AflATP 


-20 


P-E1.X 





fJ.ElX 


-5 


pEl.Xt 





PA.X* 





pAX* 


-5 


PA.X 





PE2.X 





PE2X 


-5 


* See Appendix A 



TABLE I: Simulation parameters 




FIG. 1: A PdPC with a single kinase enzyme shows bimodal 
distribution of the substrate, (a) Illustration of the PdPC. 
The kinase molecule has two slowly converting conformers 
with different enzyme activity, (b) The bulk response curve 
shows sigmoidal and monotonic zeroth-order ultrasensitivity. 
These results are obtained by solving the rate equations with 
the parameters given in Table 1. (c) A trajectory of the num- 
ber of X with a single two state kinase. The total substrate 
number Nxt ~ 100, conversion rate constant between two 
kinase states a — 0.001, and other parameters are shown in 
Table 1. (d) The distribution of X corresponding to (c) shows 
bimodal distribution. 



B. Multi-enzyme systems 

For N independent two state kinases the probability 
to have k kinases in conformer 1 has following binomial 
distribution 

7r(fc,t)= Q)p^(l-pi)^-', (6) 

where we defined 

p^{t)=pl-{pl-p[)e-^'^+P^\ (7) 

where p\ is initial probability to be in the state 1, and the 
steady state probability pf is given by pi = P/{a + (3). 

Note that binomial distribution Eq. © becomes nor- 
mal in the case A*" — > 00. Therefore, in this limit one 
can represent overall concentration of the enzymes as (cf. 
with equation S8 in [29j ) 

Et^E + m, (8) 

where ^(t) is Gaussian noise with correlation function de- 
caying exponentially fast, see Eq. ([7]). Hence, only when 
the switching rates a and /3 are fast {i.e., no dynamic 
disorder) the white noise approximation used in refer- 
ence [2^ is expected to work well, since in this case ex- 
ponential decay of the correlation function can be safely 



(theory) 



50 100 



(simulation) 



50 100 



d) 0.02 




b) 

i= 0.04 

I 

g0.02 





-X 10" 



^2 
O 



0.08 





FIG. 2: Dependence of substrate distribution of a PdPC with 
a single two-state kinase enzyme on the kinase conformer con- 
version rates. Total substrate molecule Nxt = 100 unless 
specified otherwise. Simulation (s): lower, Theory (t): up- 
per. For the parameters chosen, the theoretically predicted 
two critical points where the two peaks (singular points in 
the theory) disappear (Eq. ([5])) are ai = 0.064, ci2 ~ 0.137. 
(a) at = as — 0.001 << qi,Q2. (b) at — 0.055 close to 
Qfi. Corresponding as — 0.03. (c) at — 0.064 = ai. Cor- 
responding as = 0.04. (d) at = 0.137 — a^. Corresponding 
Os — 0.1. (e) at = as — 0.5 >> ai,a2. (f) The num- 
ber of peaks v.s. a (red stars). The vertical lines indicate 
the critical a values of peak disappearance from simulation 
(blue solid line) and theoretical prediction (red dashed line). 
The left blue solid line represents the disappearing of the first 
peak and the right blue solid line represents the disappear- 
ing of the second peak, leaving the distribution a single peak 
distribution, (g) The simulated (stars with dashed line) and 
theoretically predicted (solid line) dependence of the critical 
a value on the total substrate number Nxt- Numbers of all 
other species and the volume are increased proportionally to 
keep all the concentrations constant. 
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FIG. 3: A PdPC with 50 two-state kinase (and phosphatase) 
enzymes and 1500 substrates. Other parameters are the same 
as in 1-enzyme case, the a values are: dotted line 4 x 10~^, 
dash-dot line 0.0004, dashed line 0.02, and solid line 0.2. 
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FIG. 4: Dwelling time distribution of a PdPC with a single 
two-state kinase enzyme (q = 0.001). (a) A typical trajec- 
tory with steps indicated (dark solid line). The initial and 
final times of one dwelling state are indicated as ti and t/, 
respectively, (b) The dwelling time distribution and exponen- 
tial fitting of the upper (left panel) and lower (right panel) 



states. The fitting slopes are 
respectively. 
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gets 



{n{t)n{t')) = ^nn'P{n,t\n',t')P{n',t'), 



(9) 



replaced by 5-function. Recently Warmfiash et.al also 
discussed the legitimacy of using the (5-function approxi- 
mation (39j . 

Let us calculate noise-noise correlator explicitly. One 



where P{n, t\n', t') is the conditional probability to have 
n enzymes in conformer 1 at time t, provided that at 
time t' the number of enzymes in conformer 1 is exactly 
n' . The second term in the product, P{n',t'), is the 
probability to have n' enzymes in conformer 1 at time t' 
and it depends on initial conditions. However, for late 
times t > t' {a + (3)~^ it can be safely replaced by the 
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FIG. 5: Multi-State Enzyme fluctuation produced bimodal 
distribution, (a) Enzyme conformation fluctuation along a 
harmonic potential (blue line) and the Boltzmann distribu- 
tion (gray bars and red line), (b) The corresponding bimodal 
distribution of reactant X. 

steady state distribution: 

Pin',t')^n,{n'). (10) 

As for conditional probability, we derive 

P{n,t\n',t') ^TT{n,t~t'), (11) 

where dependence on the number n' comes from the ini- 
tial conditions in Eq. ([7]): 

Pi = (;;;). (12) 

It means that pi{t — t') in Eq. ([7|) depends imphcitly on 
n', and hence will be defined as pi{n', t — t'). After some 
algebra we get 

(n(i)n(t')) =NJ2 ^'Pii^'^ t - t') (]^) (Pi)"' (1 - Pi)"^"' 

(13) 

It is easy to check that the first term in the expres- 
sion for pi{n',t — t'), namely pf, adds the contribution 



{n{t)) {n{t')) for a late times. Therefore, one obtains a 
two time correlation function 

{n{t)n{t')) - {n{t)){n{t')) = (j2g-(a+^()(t-t')^ (^4) 

exponentially decaying in time with correlation strength 
a that can be explicitly obtained from the Eq. 

V. DISCUSSION AND CONCLUDING 
REMARKS 

Slow conformational fluctuations have been suggested 
to be general properties of proteins, and result in dy- 
namic disorder. However, so far only metabolic enzymes 
have been directly examined at the single molecule level. 
If demonstrated, the existence of dynamic disorder in 
general may greatly modify our understanding of dy- 
namics of biological networks (e.g., signal transduction 
networks). It provides a new source of in general non- 
white noises. In this work, we exploit the ultrasensitivity 
of a PdPC (or a similar system as discussed previously) 
to amplify molecular level slow conformational fluctua- 
tions. The method may be used experimentally for quick 
screening and qualitative/semi-quantitative estimation of 
molecular fluctuations in signal transduction networks. 
Here we propose a possible experimental setup. One 
adds one or a few kinases, corresponding phosphatase 
(with its amount adjusted so the average activity ratio 
6* ~ 1), a relatively large amount of substrate molecules, 
and ATP regeneration system in an isolated chamber. 
Experimentally one may consider the micro-fabrication 
technique produced high density small reaction cham- 
bers used previously in single molecule protein motor 
and enzyme studies 0, IS]. Containers may stochasti- 
cally contain different number of molecules of the kinase, 
phosphatase, with some of them giving the desired ~ 1. 
Monitoring substrate fluctuations {e.g., through fluores- 
cence) may reveal information about molecular level fluc- 
tuations. In general protein fluctuation is more compli- 
cated than the two-state model used here. The latter 
should be viewed as coarse-grained model. In this work 
for simplicity we didn't consider possible conformational 
fluctuations of the phosphatase and even the substrate 
(which may act as enzymes for other reactions). Includ- 
ing these possibilities make the analysis more difficult, 
but won't change the conclusion that molecular level fluc- 
tuations can couple to fluctuations at network levels and 
be amplified by the latter. 

There are several studies on systems showing stochas- 
ticity induced bimodal distribution without determinis- 
tic counterpart. [13, IM IH, S, S, SI In eukaryotic 
transcription, a gene may be turned on and off through 
binding and dissociation of a regulating protein, which 
may result in bimodal distribution of the expressed pro- 
tein level. The process is mathematically equivalent to 
the problem we discussed here. Physically the mecha- 
nism of generating a bimodal distribution is trivial. The 
system (PdPC) has a fluctuating parameter, the ratio 
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of the overall enzyme activity 9 (not 9) . When the pa- 
rameter fluctuates sufficiently slow, the distribution is 
approximately a mixture of localized distributions cor- 
responding to different parameter values, and thus may 
have more than one peak. This situation is fundamen- 
tally different from macroscopic bistable systems, which 
have more than one steady state for a given set of param- 
eters, and usually some feedback mechanism is involved. 
Possible biological significances of a network generating 
bimodal distributions without deterministic counterpart 
has been suggested in the literature [1^, HI], \^ . It re- 
mains to be examined whether the mechanism discussed 
in this work is biologically relevant, or reversely evolu- 
tion has selected signal transduction proteins showing 
minimal dynamic disorder jsBj • As shown in [2|; and [3l 
with fixed enzyme molecular property, a system reduces 
to unimodal distributions on increasing the system size. 
Therefore the mechanism of dynamic disorder induced bi- 
modal distribution plays a significant role only for small 
sized systems. We want to point out that Morishita et.al. 
(4^ has theoretically suggested that signal transduction 
cascades have optimal performance with only ~ 50 copies 
per specie, which makes the dynamic disorder mechanism 
plausible. In a real system it is more likely that noises 
arising from dynamic disorder, which has broad time- 
scale distribution, will couple with sources from other 
processes, such as enzyme synthesis and degradation, 
and may result in complex dynamic behaviors. There- 
fore physical chemistry studies of molecular level protein 
dynamics may provide important and necessary informa- 
tion for understanding cellular level dynamics. 



APPENDIX A: APPENDIXES 



Here we show how we make connections between the 
equations for bulk analysis and the ones for stochastic 
simulations with molecular number. For example, if we 
have a reaction 



X 



El ^ X El . 

"•1 (, 



we can write down the ODE equations for this reaction 

First of all, we choose as our time unit, where k-f 

is the rate constant for the backward enzymatic step. 
Then, 

'^^^^ = k,f[x]m-k,,[xE{\. 

with kif = k^j/k-f, kih — fcjj/fc_j. If we want to deal 
with variables in the unit of molecular numbers instead 



of concentration, we then have 



dt 



k 



Nx N. 



1/ 



N 



16- 



NaVoNaVo '"NaVo 



where Na is Avogadro constant. Vo is the volume of the 
system. We can further simplify the expression, 

dNxE, , NxNe, , 



dt 



NaVo 



In all of our simulations, we kept a constant value for the 
substrate concentration Nxt/iNAVo) = f • Then, 



dNxEi _ hf 



dt 



N 



NxNe, -kuN_ 



Ib-l^XEi 



Xt 



B: Multi-state enzyme fluctuation 

In general, an enzyme fluctuates continuously along 
conformational coordinates. One should consider the 
two-state model discussed in the main text as a coarse- 
grained model. Here we will use a more complicated 
model to show that our main conclusions still hold in gen- 
eral case. We consider an enzyme diffuse slowly along a 
harmonic potential of coordinate x, G{x) = x^, where we 
have chose the units so G = IksT at a; = I. Motion along 
the conformational coordinate couples to the enzymatic 
reaction rate constants with an exponential factor fc+y = 

^ exp(Aa;), where A = —0.5. One specific example is 
that X is the donor-acceptor distance for an electron- 
transfer reaction. We model the diffusion as hopping 
among 10 discrete states. The conversion rate constant 
between enzyme conformations is at 10~^. The back- 
ward rate constant is then determined by the detailed 
balance requirement Pi = aiexp(— (xf — xf_^)/kBT). 
Other parameters are the same as in the 2-state enzyme 
simulations. The reactions are listed below. 



X -f- El ^ X El ^ X 

kib *:+!(> 

X + E2 ^ X E2 ^ X ' 

k2b k+2b 



El, 
E2, 



fe(m-l)/ fc+(m-l)/ 

X + Em-1 ^ XE„i-i ^ X +£",„_!, 

km/ fe+m/ 

X + Em ^ XEjn ^ X* + E„i, 
kmb k^,„b 

kaf k^ f 

X* + A ^ X*A ^ X + A, 

kab k^b 

ai 0.2 OLm-2 
El ^ E2 ^ ^ Ejn-1 ^ Em, 

Pi th Pm-2 Pm-l 

EiX ^ E2X ^ ^ Em-lX ^ EmX. 

01 P2 PL-2 P'm-1 

O shows that the reactant X shows a bimodal distribu- 
tion if the enzymatic conformational fluctuation is slow. 
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This result reiterates our suggestion that the ultrasensi- 
tive network amplifies small enzymatic activity fluctua- 
tions into large substrate number fluctuations. 



By omitting the diffusion terms in the Eq. ([2]) and Eq. 
([3]) one derives for steady state: 



dx [fi{x)pi\ = api - [3p2, 
fi{x)pi + f2{x)p2 = const. 



(Al) 
(A2) 



For very fast switching rates a and (5 we expect an uni- 
modal distribution centered somewhere in between two 
'turning' points xi and X2, see Eq. Q. Therefore, the 
steady state solution should vanish at these points and 
be identically zero outside an interval [0:2,2:1]. 

In order to satisfy these boundary conditions, one has 
to set a constant in Eq. (|A2p to be zero. Hence, we 
obtain 



dx [h{x)pi] = 



h{x) 



Pi- 



rn 



The solution of the Eq. jASj) depends, of course, on par- 
ticular choice of the function f{x). However, it is guar- 
anteed that there is an unique root of the function fix ) 
in the corresponding physical region of variable x [29[. 
Hence, the differential equation Eq. (|A3p is singular only 



at two points X2 and Xi, which are the boundary points. 
One can find an asymptotic behavior of the steady state 
solution near these points: 



Pi - {x-X2THxi-x)'''g{x), 



(A4) 



where g{x) is analytic function of a; in the interval [x2, Xi], 
satisfying condition g{x2) = g{xi) = 1. The exponents 
02 and ai are 



02 



ai 



f2{x2) 

a 



-1, 
- 1. 



(A5) 
(A6) 



Therefore, if the conditions Eq. ([5]) are satisfied, one 
expects unimodal distribution. Otherwise, there exists at 
least one additional peak in the distribution. In this case 
of the slow switching rates the equation Eq. (|A3p predicts 
divergence of the solution at one or both boundary points 
X2 and Xi- This is an indication that diffusion terms 
that we omitted for our estimate become relevant. The 
diffusion terms makes the overall distribution finite. 
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